Mon. Not. R. Astron. Soc. 000, 000-000 (0000) Printed 2 May 201 1 (MN KTeX style file v2.2) 



A comparison of spectroscopic methods for detecting starlight 
scattered by transiting hot Jupiters, with application to Subaru data 
for HD 209458b and HD 189733b 

(N Sally V. Langford, 1,2 J. Stuart B. Wyithe, 1 Edwin L. Turner, 2 3 Edward B. Jenkins, 2 
b^; Norio Narita, 4 Xin Liu, 2 Yasushi Suto 2 ' 5,6 and Tom Yamada 7 

1 School of Physics, University of Melbourne, Parkville, Victoria 3010, Australia 



Cs 
(N 

Oh 
i 

o 



(N 
> 
(N 
OS 

m 

SO 
O 
O 



2 Princeton University Observatory, Princeton NJ 08540, USA 

3 Institute for the Physics and Mathematics of the Universe, University of Tokyo, Kashiwa 277-8568, Japan 
^National Astronomical Observatory of Japan, 2-21-1 Osawa, Mitaka, Tokyo, 181-8588, Japan 
5 Department of Physics, The University of Tokyo, Tokyo 113-0033, Japan 

^Research Center for the Early Universe, School of Science, University of Tokyo, Tokyo 113-0033, Japan 
' Astronomical Institute, Tohoku University, Sendai 980-8578, Japan 



2 May 2011 



ABSTRACT 

The measurement of the light scattered from extrasolar planets informs atmospheric and 
formation models. With the discovery of many hot Jupiter planets orbiting nearby stars, this 
motivates the development of robust methods of characterisation from follow up observations. 
In this paper we discuss two methods for determining the planetary albedo in transiting sys- 
tems. First, the most widely used method for measuring the light scattered by hot Jupiters 
dCollier Cameron et alJ [T999) is investigated for application for typical echelle spectra of a 
transiting planet system, showing that a detection requires high signal-to-noise ratio data of 
bright planets. Secondly a new Fourier analysis method is also presented, which is model- 
independent and utilises the benefits of the reduced number of unknown parameters in tran- 
siting systems. This approach involves solving for the planet and stellar spectra in Fourier 
space by least-squares. The sensitivities of the methods are determined via Monte Carlo sim- 
ulations for a range of planet-to-star fluxes. We find the Fourier analysis method to be better 
sui ted to the ideal case of typica l observations of a well constrained transiting system than 
the ICollier Cameron et al.l 0999) method. To guide future observations of transiting planets 
with ground-based capabilities, the expected sensitivity to the planet-to-star flux fraction is 
quantified as a function of signal-to-noise ratio and wavelength range. We apply the Fourier 
analysis method for extracting the light scattered by transiting hot Jupiters from high resolu- 
tion spectra to echelle spectra of HD 209458 and HD 189733. Unfortunately we are unable to 
improve on the previous upper limit of the planet-to-star flux for HD 209458b set by space- 
based observations. A lcr upper limit on the planet-to-star flux of HD 189733b is measured 
in the wavelength range of 558.83-599.56 nm yielding e < 4.5 x 10~ 4 . This limit is not 
sufficiently strong to constrain models. Improvement in the measurement of the upper limit 
of the planet-to-star flux of this system, with ground-based capabilities, requires data with a 
higher signal-to-noise ratio, and increased stability of the telescope. 

1 INTRODUCTION with an expected flux of less than 10~ 4 times the direc t stellar flux 

dCollier Cameron et al.lll999l ; ICharbonneau et al.lll999h . 

With increasing numbers of hot Jupiter planets currently being 
More than 500 planets have been discovered outside the Solar Sys- discovered, follow up observations become very important for char- 
tern in the past two decades. Among these planets is a class of ob- acterisation that requires longer observation times, and for measur- 
jects known as hot Jupiters. These are rough ly Jupiter-mass ob - ing the accuracy of the orbital parameters. The space-based CoRoT 
jects that orbit their host stars within 0.05 AU ( Seager et alJuOOOf) . and Kepler missions, which are currently monitoring many thou- 
but not close enough to be considered very hot Jupiters, where sands of nearby stars, will yield a wealth of new transiting planets, 
thermal emission dominates the visible scattered light. It has been These planets will also be followed up with radial velocity mea- 
predicted that hot Jupiter planets will allow the starlight scattered surements to constrain the orbital parameters. In order to fully re- 
by their upper atmospheres to be directly detected in the optical, alise the scientific potential of these systems, robust methods of 



2 Langford et al. 



analysis of the starlight scattered by the transiting hot Jupiter plan- 
ets are required, as well as an understanding of their application to 
various systems. 

There has been no definite detection of the broadband un- 
polarised visible light scattered by a hot Jupiter. Extrasolar plan- 
ets have only recently been directly imaged in the optical. Th ese 
planets orbit their parent stars at ~100 AU dKalas et al]|2008h . In 
addition, visible polarised scattered light has been detected from 
the extended atmosphere of the transitin g planet HP 189 733b 
( iBerdvugina et alj2008|]201ll) . However. Iwiktorowiczl d2009h was 
unable to reproduce this result. Measuring the light scattered by 
a hot Jupiter planet is challenging due to the small angular sepa- 
ration, and the small planet-to-star flux ratio, which can be of the 
order of the systematic noise of ground-based observations. Ex- 
tracting the planet signal from the combined noisy stellar and scat- 
tered flux therefore generally requires making assumptions about 
the system. Upper limits on the visible light scattered by a hot 
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l20ld : ISing & L6pez-Moralesll2009h . 

The light scattered by an extrasolar planet will depend on the 
composition of the atmosphere as well as the orbital parameters , 
allowing the test ing of atmospheric models dMarlev etal] 1 19991 : 
iGreen et al.|[2003b . Models predict a range of planet-to-star fluxes 
for different classes of planetary atmospheres, which can guide the 
search f or thermal emission and re flected light from hot Jupiter 
planets dSudarskv et ai]|200d 120031) . The non-detections of scat- 
tered light and resulting upper limits on the albedos of hot Jupiter 
planets aids the development of theories of extrasolar planet at- 
mospheres. The lack of a detection of starlight scattered by the 
planet HD 209458b suggests that the atmosphere is much less 
reflecti ve than previously expected, and lacks highly scatterin g 
clouds dSeager et al.ll2000l:lRowe et alJl2006l ; lBurrows et alj|2008h . 
Increasingly sophisticated models are attempting to explain the in- 
flated radius and absence of a verifiable detection of light scattere d 
by this hot Jupiter planet dHood et alj|2008l ; iBurrows et ai]|2008T) . 
Further investigation and observation of known hot Jupiter planets 
will continue to guide the models, and potentially yield a definite 
detection of scattered broadband unpolarised visible light. 



1.1 Current Methods 

Prior to the discovery of transiting extrasolar planets, two meth- 
ods of measuring the planet-to-s tar flux of hot Jup i ter p lanets 
were concurrently develope d by ICharbormeau et alj (l999) and 
Collier Cameron et a i] dl999h . Both methods require modelling the 
direct stellar contribution, which is then removed, leaving the 
planet signal buried in noise. They model the system over a range of 
orbital phases, inclinations and velocities in order to find the signa- 
ture of the planet's scattered light, and the likelihood of the best fit 
orbital parameters. Both methods been imp ortant in beginning to 
const rain models of planetary atmospheres ( Sudarsk v et al.ll2000l 
120031) . 

An upper limit on the planet-to-star flux of the extrasolar 
planet r Bootis b was measured to be < 5 x 1 times the direct 
stella r flux at the 99 per cent confidence level dCharbonneau et al.l 
1999). This measurement corresponds to a wavelength range of 
465.8^-98.7 nm, and assumes an inclination of i > 70°. The result 
was limited by the photon nois e of the data, not systematic e ffects. 

In disagreement with this, ICollier Cameron et alj < fl999h mea- 



sured a best fit to the planet-to-star flux of the same system to be 
(7.5 ± 3) x 10" 5 at the 97.8 per cent confidence level. This de- 
tection contains systematic uncertainty due to some dependence on 
wavelength of the albedo of r Bootis b. For the wavelength re- 
gion of 456-524 nm the planet-to-star flux ratio was found to be 
(1 .9 ± 0.4) x 10~ 4 , which is brighter than the upper limit set 
bv lCharbonneau et alj(fl999h . Outside this wavelength region there 
was no signal detected. The wavelength dependence of detected 
scattered light is useful in investigating the classes of hot Jupiter 
planet atmospheres. Disagreement in planet-to-star flux measure- 
ments may be due to the different phase functions adopted for the 
planets, and the methods used for subtracting t he stellar contribu- 
tion. For example. lCollier Cameron et al. I d 19991) adopt a Venus-like 
phase function of the planet to account for higher backscattering 
from clouds, as compa red to a simple Lambert-law phase func- 
tion. lLeighetal.lf2003ah failed to confirm the tentative detection of 
light scatter ed by r Bootis b with subseq uent reanalysis and addi- 
tional data. ICollier Cameron et alj d2000h later revised their result, 
setting an upper lim it on the planet-to-star flux of<3.5xl0~ . 
iRodler et al] j201Ch recently measured an upper limit on the light 
reflected by r Bootis b to be < 5.7 x 10 -5 times the direct stellar 
flux, at the 99.9 per cent signifi cance level. This upper l imit was 
based on the method devised by 
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ICollier Cameron et alj d2002h set an upper limit on the v An- 
dromeda system with a 0.1 per cent false alarm probability of 
< 5.84 x 10~ 5 . With the same method, iLeigh et alj d2003bl) 
placed an upper limit on the planet-to-star flux of H D 75289b of 
<. 4.1 8 x 10~ 5 at the 99.9 per cent level. However. IRodler et all 
(2008) disputed this result, stating that it was based on incorrect 
orbital phases for the planet. They rederived the upper limit on the 
planet-to-star flux of HD 75289 b to be ~ 60 p er cent higher, via 
the approach taken by Char bonneau et al] dl999l) . 

These methods yield varying results for the same systems, 
most likely due to the difficulty in measuring the dim planet-to- 
star flux, their model dependence and the unconstrained orbital pa- 
rameters. However, the up per limits set are lower than p redicted by 
theoretical considerations dCollier Camer on et al . 1999). This sug- 
gests that hot Jupiter planet atmospheres are less reflective than ex- 
pected based on the solar system gaseous planets. Constraining at- 
mospheric models further requires deeper upper limits on the scat- 
tered light, and minimising the model-dependence of measuring the 
albedos of hot Jupiter planets. 

I n order to r educe the parameter space and number of assump- 
tions, iLiu et al. I d2007l) developed a method specific to detecting 
the starlight scattered by transiting hot Jupiter planets. Knowing 
the planet's velocity from orbital constraints allows a shift-and-add 
method to be adopted, to verify the presence of the scattered light. 
Using the equivalent wi dth ratios of the scattered and direct com- 
ponents of the spectra, ILiu et al. I d2007l) found the planet to star 
ratio of HD 209458b to be (1.4 ± 2.9) x 10" 4 in the wavelength 
range 544-681 nm. This result is expected to be updat ed, following 
a corr ection for the non-linearity of the HDS CCDs dTaiitsu et al.1 

boich . 

Space-based satellites can use transit photometry to mea- 
sure the starlight scattered by hot Jupiters at a higher signal-to- 
noise ratio than possible with ground-based capabilities. To de- 
termine the amount of starlight sc attered by HD 209 458b from 
the depth of the secondary eclipse, iRowe et all d2008l) measured 
light curves of HD 209458 for the duration of the planetary orbit 
with the Microvariability and Oscillations of Stars (MOST) satel- 
lite. Their observations yielded a la upper limit on the planet- 
to-star flux ratio of < 1.6 x 10 -5 . In the wavelength region of 
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400-700 nm this corresponds to a geometric albedo of 0.038 ± 
0.045 which is much less reflective than the solar system's gi- 
ant planets. This rules out bright clouds at a high altitude, and 
is consistent wi th the low albedo upper limits of v Andromed a 
and HD 75289 dCollier Cameron et all [20021 : llxigh et alj|2003bl) . 
While searching for transiting planets, the CoRoT and Kepler 
missions have also yielded a range of planetary upper limits on 
the albedos of the hot Jupiter planets su ch as Kepler-7b, Kepler- 



5b, CoRoT-lb, CoR o T-2b and OgleTr56 ( |Kipping & B akos 
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2011 



2010; 



Sing & Lopez-Mo rales 2009]). Further studies and repeated obser- 



vations of these bright, transiting systems will continue to constrain 
atmospheric models and theories of the formation of solar systems. 
However, space based missions are costly and it is not feasible to 
consistently follow up on the expected yield of hot Jupiter plan- 
ets while the search continues for Earth-like extrasolar planets. 
It is therefore beneficial, and timely to develop a suitable model- 
independent method of extracting the planet-to-star flux from tran- 
siting hot Jupiters in ground-based observations. It is also advan- 
tageous to measure phase independent planet-to-star fluxes with- 
out assuming a grey albedo. This motivates a model-independent 
Fourier based method for measuring the light scattered by hot 
Jupiter planets, which we introduce in this paper. This method re- 
quires minimal assumptions about the system and is suited to tran- 
siting planets where parameters are constrained by survey observa- 
tions. 

In this paper we quantify th e precision of the existing method 
of ICollier Cameron et al. (1999) in measuring the planet-to-star 
flux of transiting hot Jupiter planets via Monte Carlo simulations. 
We also introduce and quantify a new method specifically tailored 
for application to transiting planets in order to constrain the orbit 
and detect smaller signals with ground-based capabilities. We begin 
with a review of the method o f measuring the light scattered by ex- 
trasolar planets developed by ICollier Cameron et alj (l999). This 
method adopts a detailed approach to extracting the signal of the 
planet from the residual noise of the spectrum by least-squares de- 
convolution, following the stellar spectrum removal. Monte Carlo 
simulations are adopted to determine the sensitivity of this method 
for application to mock high resolution echelle spectra of a transit- 
ing hot Jupiter planet system, observed with ground-based capabil- 
ities. A new method is then developed which is model-independent 
and utilises the benefits of reduced parameters for a transiting sys- 
tem. This Fourier analysis method involves solving for the planet 
and stellar spectra in Fourier space by least-squares, and is applica- 
ble to transiting hot Jupiter planets. To guide future observations, 
the expected sensitivity to the planet-to-star flux fraction is quan- 
tified as a function of signal-to-noise ratio and wavelength range. 
This method is then tested on Subaru HDS spectra of HD 209458 
and HD 189733, which host a transiting hot Jupiter planet. 



2 SYSTEM PARAMETERS 

A planet will reflect an amount fp of its host star's flux f, 

2 



W^A) fRp(A) 
e(a, A) = r = p(A)g(a, A) 1 



(1) 



where a is the semi-major axis, Rp(A) is the radius of the planet, 
p(A) is the geometric albedo, and g(a, A) is the phase func- 
tion a t planetary phase angle a and wavelength A jLeigh et"al] 
l2003ch . From theoretical considerations a planet-to-star flux of 



e ~ 10 4 is expected fo r a hot Jupiter planet at full phase 
dCollier Cameron et al.ll 19991) . 



2.1 Phase Function 

The flux scattered by a hot Jupiter planet is a phase dependent frac- 
tion of the maximum reflectance at full phase. Unfortunately, due 
to a lack of signal and phase coverage it is not yet possible to de- 
termine the exact phase function g(cv) for a hot Jupiter planet being 
observed, and if required, it must therefore be modelled as a func- 
tion of the phase angle a. 

It is simplest to adopt a Lambert-law sphere which assumes 
isotropic scattering over 2-7T steradians. This typically has the form; 



»(a) = [sin a + (n — a) cos a]/n, 



(2) 



where the phase angle a is determined from the inclination i and 
the orbital phase (\>\ 



cos a = - sin i cos 2n<f). 



(3) 



For an orbital phase of (f> — the planet is closest to the observer, 
and has the smallest planetary phase. For an orbital phase of (f> — 
0.5, the planet is at the furthest point in its orbit from the observer, 
with maximum reflectance. 

To account for stronger back-scattering than the Lambert-law 
phase function, it may be more realistic to adopt the phase function 
of a similar planet in the Solar System, that has a cloud-covered 
surface. For example; 



g(a) = 10 
where; 



-0.4Am(a) 
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0.09 
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+ 2.39 



100° 



0.65 
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(4) 



(5) 



which is a polynomial appr oximation to the empirically determined 
Venus-like phase function dCollier Cameron et al]|2002l) . For tran- 
siting planets, which have orbital inclinations close to 90°, these 
two functions have a similar variation in flux of less than 10 per cent 
over the orbital phases adopted for the mock spectra presented here, 
and the Subaru HDS spectra of HD 189733b and HD 209458b. 



2.2 Velocity 

A planet's velocity relative to its host star Vp (<f>) at an orbital phase 
<f> is given by; 



Vp(</>) = Kpsin(27r</>), 



(6) 



where tf> = (t — To)/P or tj at the time t, where To is the transit 
epoch and Pqj-]-, is the orbital period, for a circular orbit (ellipticity 
e=0). The apparent radial-velocity amplitude Kp about the centre 
of mass of the system is; 



Kn 



27ra sini 
Porbl + q' 



(7) 



where a is the orbital distance, i is the inclination and q = Mp /M* 
is the mass ratio of the planet to the star. For a system under con- 
sideration, a, 2 and Rp can be determined from radial-velocity mea- 
surements and transit photometry. For non-transiting planets these 
parameters have greater uncertainty, and Kp is determined from fit- 
ting observational light curves. 
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Figure 1. Mock spectra constructed from high signal-to-noise ratio solar spectrum. Panel A - Example mock spectra for red wavelength 
range of 560-608 nm. Panel B - Example mock spectra for blue wavelength range of 440-515 nm. Statistical Gaussian noise added to the 
spectrum is distributed with standard deviation \/N for N counts per pixel. The wavelength is logarithmically binned such that the Doppler 
shifted scattered light can be added as an integer pixel shifted copy of the star, with amplitude scaled by eg(a). The orders of the spectra were 
combined and normalised with an 1 1-knot cubic-spline. 



3 CONSTRUCTION OF MOCK SPECTRA 

Mock spectra were constructed in order to quantify the precision 
of methods for measuring the light scattered b y extrasolar planets, 
including that bv lCollier Cameron et alj d 19991) . as well as the new 
Fourier analysis method, introduced in Section 5. Spectra were con- 
structed from high resolution solar spectra (A/ A A ~ 300,000), 
obtained with the Fourier Transform Spectrometer at the Mc- 
Math/Pierce Solar Telescope situated on Kitt Peak, ArizoncQ. The 
mock spectra were constructed to be analogous to observing a 
bright hot Jupiter transiting system, such as HD 209458, with typ- 
ical ground-based 8 m telescope capabilities, such as Subaru HDS. 
An example mock spectrum is shown in Fig.[TJ The high resolution 
solar spectra were scaled to the counts expected from HD 209458 
with Subaru HDS for a typical exposure time of 500 seconds. A 
spectral resolution of R = A/AA = 45,000 was mimicked by 
Gaussian smoothing, consistent with twice the critical oversam- 
pling. Noise due to statistical fluctuation in photons from the hy- 
pothetical measurement with Subaru HDS was added after includ- 
ing the planet contribution and smoothing. The noise is Gaussian 
distributed with a standard deviation of \/N for N counts per pixel. 
The planet velocity dictates the wavelength shift of the scattered 
spectrum, taken as a scaled copy of the stationary stellar spectrum. 
During the length of an exposure the planet velocity varies by less 
than the width of the absorption features. 

Planet velocities were randomly selected within the range 30- 
80 km s _1 towards the observer. This is consistent with the range 
of velocities of HD 209458b when the planet is just out of sec- 
ondary eclipse and the orbital phase is in the range 0.55-0.60. 
There is less than 10 per cent variation in the flux scattered from 
the planet within this phase range based on either a Lambert- sphere 
or a Venus-like phase function. This is a benefit of using transiting 
systems to measure planet-to-star fluxes. As the planet velocity am- 
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plitude can be measured from the constrained parameters, a smaller 
range of the orbit can be observed and selected to be where the 
planet is close to full phase and maximum reflectance. The veloc- 
ities must be large enough to shift the planet signal clear of the 
stellar absorption lines. 

Instrumental variations, such as flexure of the spectrograph 
and guiding errors, can cause wavele ngth variations be tween the 
spectra during a night of observations (Winn et al. 2004). This can 
be accounted for in data reduction by adjusting the wavelength cal- 
ibration using the star's spectral absorption lines. Variation in the 
flux calibration between exposures is removed by continuum fitting 
and normalisation. The mock spectra are taken to be ideal, with no 
variation in the shape of the response of the CCD with time and 
the host star spectra in all exposures are assumed to be perfectly 
aligned. 



4 COLLIER CAMERON METHOD 

The most widely used current method ( henceforth the Collier 
Camer on method) was first pr esented by ICollier Cameron et al.l 
1 19991), and further developed bv lCollier Cameron et al. 1 2002h and 
iLeigh et all (2003b). The Collier Cameron method has previously 
been adopted to attempt to measure the li ght scattered by three 
different, non-transiting hot Jup iter planets I Collier Cameron et al.l 
1 19991 12OO2I ; ILeigh et alj[2003brt . It involves detailed modelling of 
the stellar contribution and matched-filter analysis to extract the 
planet signal from the noise. It is a useful method for the case 
where the hot Jupiter planet's orbit is unknown, and has allowed 
constraints to be placed on highly reflective planetary atmospheres 
dRodler et alj|2010h . In this Section, the Collier Cameron method 
is reviewed for application to a transiting planet system, using the 
mock spectra described in Section [3] within Monte Carlo simula- 
tions. 

For the analysis of non-transiting planets, the Collier Cameron 
method requires observations of as much of a full orbit as possible, 
such that the most likely orbital phase function and planet velocity 
can be measured. For transiting systems the velocity of the planet 
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Figure 2. Example probability distribution. An example probability dis- 
tribution with a mean planet-to-star flux of e = 1.2 X 10 -4 marked by the 
vertical solid line. The distance of one, two and three standard deviations 
from the mean are marked by the vertical dashed lines. This distribution has 
a statistical significance of > 1<t as the mean is greater than one standard 
deviation, but less than two standard deviations from zero. 

is known, and a small section of the orbit can be observed. As we 
are evaluating the Collier Cameron method for the case where the 
planet is transiting, the thirty mock spectra that were used to mea- 
sure a planet-to-star flux are in a small orbital phase range, either 
side of the secondary eclipse, when the majority of the dayside 
of the planet is facing the observer. This corresponds to velocities 
with magnitudes of 30-80 km s^ 1 . The planet's spectral features 
are therefore shifted clear of the corresponding stellar features, and 
the planet is close to maximum brightness. The planet is assumed 
to follow a Venus-like phase function as defined in Section [2~T1 
The initial input spectra do not need to be merged for this method, 
but can remain as separate orders, without continuum normalisa- 
tion. The spectra are in the red wavelength range of 560-680 nm, 
and the blue wavelength of range of 440-515 nm. We note that the 
mock spectra are ideal, in that they are all aligned, without cosmic 
ray hits or vignetting, and the stellar contribution is constant. 

The Collier Cameron method requires modelling and remov- 
ing the direct starlight from the total spectrum, leaving the planet 
signal buried in noise. It is very difficult to accurately remove the 
stellar contribution, as small distortions in the strong stellar lines 
can produce changes larger than the planet signature. Problems can 
arise due to shifting of the spectrum on the detector, changes in 
the telescope focus and atmospheric seeing. In the blue wavelength 
range, the density of strong lines makes removing the stellar sig- 
nal particularly difficult, and observations for previous upper limit 
measurements have generally been in a central wavelength range. 

In this paper, results are presented for the ideal mock spectra, 
and the stellar spectrum does not vary with exposure. This aids the 
removal of the modelled stellar spectrum, since we do not have to 
account for sources of distortion. This simplification of the removal 
of the direct light is equivalent to perfectly modelling the stellar 
spectrum, and is th e best case scenario for the Collier Cameron 
method described in ICollier Cameron et al.1 f20o3) . The sensitivity 
of the method to the planet-to-star flux will be impeded by instru- 
mental variations, and the difficulty in modelling the stellar spec- 
trum, particularly in the blue wavelength. 

The planet signal strength is extracted from the residual after 
the direct flux subtraction via least-squares deconvolution with a 
set of known weig hts of stellar absorption lines. The method is de- 
scribed in detail in lCollier Cameron et alj |2002). This results in a 
deconvolved velocity profile which represents the lines present in 
the planet spectrum. It is similar to a cross-correlation, but neigh- 
bouring blended lines are dealt with more accurately. For the mock 
spectra, a library of spectral lines in the required wavelength range 
was compiled from d ata for the Sun fr om the Vienna Atomic Line 
Data Base (VALD) jKupka et al.lll999h . 



There are typically ripples in the deconvolved stellar profile 
at low velocities tha t may cause spurious effects on measuring the 
planet contribution ( Collie r~Cameron et alj|2002l) . For the Monte 
Carlo simulations the planet signal was extracted from the range 
of velocities with magnitudes greater than 30 km s _1 . To extract 
the planet signal, a matched-filter is constructed, and fit to the 
data to determine the planet-to-star flux. The velocity profile of the 
planet signal is modelled as a moving Gaussian, based on the form 
used in Collier Cameron et ID §002). The planet-to-star flux is de- 
termined by scaling the matched filter to match the deconvolved 
planet line profile, and findin g the best fit with a x 2 minimisation 
(ICollier Cameron et al.lll999b . 

4.1 Monte Carlo Analysis 

The sensitivity of the Collier Cameron method was determined for 
a range of signal-to-noise ratio strengths and planet-to-star fluxes 
via Monte Carlo simulations. The typical signal-to-noise ratio per 
pixel of the mock spectra was varied by increasing the counts per 
pixel, corresponding to a longer exposure time. A larger telescope 
collecting area or observing a brighter star would also increase the 
counts per pixel. The number of standard deviations separating a 
detection from zero was determined by the mean of the distribution 
of all Monte Carlo simulation solutions, divided by the standard 
deviation. This defines the distance of the mean value from zero in 
units of the spread of the distribution. For example, a result quoted 
at la has a distribution with a mean value that is one standard de- 
viation away from zero, as shown in Fig. [2] 

In the red wavelength range of the mock spectra of 560- 
680 nm, with 1465 absorption lines, a signal-to-noise ratio of 
~ 800 is required for a la detection of a planet-to-star flux of 
1.6 x 10 -4 . The blue wavelength range of 440-515 nm has many 
more absorption lines (2570), such that the method is more sensi- 
tive to smaller planet-to-star fluxes for data with the same signal- 
to-noise ratio, assuming the stellar spectra is perfectly removed. 
Planet-to-star fluxes greater than 10~ 4 are detectable at the la 
level. A planet-to-star flux of 1.6 x 10 requires a signal-to-noise 
of around ~ 800 to be detected at the 3a level, for the ideal case. 
The results of the Monte Carlo simulations with varying signal-to- 
noise ratio for the red wavelength range of 560-680 nm and blue 
wavelength range of 440-515 nm are presented in Fig. [3] There is 
a limit on the possible signal-to-noise ratio of the spectra due to the 
capability of the CCD and variations in the planet velocity during 
longer exposure times smearing the scattered signal. 

The results of the Monte Carlo simulations for ideal spec- 
tra are consistent with detections via this method for real spec- 
tra. Probable detections have been made for planet-to-star fluxes 
brighter than e = 10~ 4 , and up per limits are set with low fa lse- 
alarm probabilities. For example. [Collier Cameron et al.1 j 19991) de- 
tected a best-fit to the signal of the light scattered by r Bootis b of 
e = 1.9 x 10~ 4 for data with a typical signal-to-noise ratio of 
~1000, in a central wavelength range. 

As the Collier Cameron method was developed prior to the 
discovery of transiting hot Jupiter planets, it did not take advantage 
of the minimisation of unknown parameters in these systems. This 
method can benefit from the known planet velocity in extracting the 
scattered signal. However in constructing the template spectrum, 
the spectra are summed to smear out the planet signal. This requires 
a large range of velocities, or isolating spectra with the planet at 
minimal illumination. Observing a small section of the transiting 
planet orbit may therefore affect the accuracy of the stellar model, 
and reduce the planet signal in the residual spectrum. 
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Figure 3. Signal-to-noise ratio for mock spectra with Collier Cameron method. Contour plot of la, 2a and 3a levels of detection (solid) 
for signal-to-noise ratio values in the range 10 2 -10 3 , and e values in the range 10 _5 -2 X 10 -4 . Panel A - Results for the red wavelength 
range 560-680 nm. Panel B - Results for the blue wavelength range 440-515 nm. The dots represent the sampling of the parameter space. 



W hile the analysis methods develop ed bv lCharbonneau et al.l 
dl999l) and ICollier Came ron et a are required for a non- 

transiting system, they have not yielded a detection of the light 
scattered by a hot Jupiter planet. With the focus of many cur- 
rent extrasolar planet surveys on transiting objects, it is therefore 
timely to develop a method of extracting the planet-to-star flux spe- 
cific to transiting hot Jupiter planets, that does not require mod- 
elling the stellar contribution. In the next Section we describe a 
model-independent Fourier analysis method that requires minimal 
assumptions about the system and is suited to transiting planets 
where the parameters are constrained by other space-based, survey 
observations. 



5 FOURIER ANALYSIS METHOD 

The Fourier analysis method for extracting the scattered signal 
from the combined planet and star flux described in this Section 
is based on a method developed to deal with the effect of the 
fixed- pattern resp onse function of a detector on a moving spec- 
trum djenkinsll2002r) . In the case of a transiting planet, the stellar 
spectrum is considered to be stationary (analogous to the response 
function), while the planet's scattered spectrum moves with each 
exposure. Separating out a moving spectrum via a Fou rier trans- 
form i s a useful technique for binary systems (see e.g. lHadraval 
( 1995)). The observed spectrum recorded by the echelle spectrom- 
eter T(A), consists of the direct starlight contribution S(A), plus the 
starlight scattered by the orbiting planet P(A). The planet spectrum 
is shifted with respect to the stationary stellar spectrum with a pre- 
dictable Doppler shift, dependent on the velocity of the planet at 
time of the exposure. There is also a contribution from statistical 
noise that is not frequency dependent. 

For the total spectrum T(A) we define: 

*.(«) = i?e{J?[T(A)]}, (8) 
ty{u) = /m{^[T(A)]}, (9) 

where J? is the Fourier transform operator. Similarly p x (u)),p y (u)), 
s x (u>) and s y (ui) are defined for the planet and stellar spectra. 
As each spectrum has a different Doppler velocity displacement 
for P(A), depending on the observed velocity of the planet, the 
zero shift complex function p x (ui) + ip y (ui) is multiplied by 



exp(— 27riaj5fc), where 8k is the magnitude of the shift of P(A) in 
the k th exposure compared to the zero velocity case. 

The minimum of the sum of the squared real and imaginary 
residuals (the minimum of Q 2 ) is solved to determine the unknown 
planet and stellar spectrum coefficients p x (ui) and p y (ui), s x (uj) 
and s v (lo). This is done by setting the partial differentials of Q 2 
defined below with respect to the four unknowns p x , p y , s x and s y 
at each uj, to be zero; 

n 

Q 2 = ^(p^cosAfc +PysinA fe + s x - t x>k ) 2 (10) 
fe=i 

n 

+ sin Afc +p y cos A fc + s y - t y , k ) 2 , (11) 

where Afc is an abbreviation for 2nu)Sk- 

Therefore, the best solutions for p x (cu), p y (uj), s x (uj) and 
s y (ui) are found by solving a system of four linear equations at 
each ui. The linear equations can be summarised as; 

C(u) ■ u(cj) = b(w), (12) 

where u(cj) is a vector of p x ,p y , s x and s y at each ui and the coef- 
ficients dj and hi are listed in TableQ] 

Singlular value decomposition is used to solve the system of 
linear equations. It can be verified that the solution is a minimum 
in Q 2 by plotting nearby solutions. After solving for each vector 
u(u>) over all ui, the inverse Fourier transforms of the two pairs 
of terms ui(ui) + iu2(ui) and 113(01) + iu4 (uj) recover the best 
representations of P(A) and S(A). The planet spectrum, P(A), is 
divided by the stellar spectrum, S(A), and the mean found over 
all wavelengths, yielding the corresponding planet-to-star flux. The 
solution can be verified to be the planet contribution by repeating 
the measurement while shuffling the order of the Doppler shifts 
with respect to their exposures, which will yield a null result. 

If the velocity displacements of the observations in a set are 
equally spaced, the the matrices become singular and the solu- 
tions can blow up. Therefore, the timing of the observations should 
be such that the velocity shifts are noncommensurate. This is not 
unique to the Fourier analysis method, and should be considered 
when planning an observing program. 

Solving for the planet and stellar spectra in this way has the 
following benefits; (i) it does not require a model of the star, (ii) 
it is still effectively using the signal of all of the absorption lines, 
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j=i 



j=2 



j=3 



j=4 



i=l n ^cosA fc — ^sinA fc t x k cos A fc — y] t yk sin A fc 

i=2 n yjsinA fc JJcosA fe ^ sin A fc + t y]fc cos A fc 

i=3 yjcosA fc JJsinAfc n X)**,* 

i=4 -yjsinA fc yjcosA fe n X)*y,* 



Table 1. Coefficients for the matrix C(uj) and vector b(u). All summations are from k = 1 to n, where n 
is the number of spectra. 



although they are not summed over to increase the signal-to-noise 
ratio, meaning that the planet-to-star flux solution can be wave- 
length dependent, and (iii) there is no assumption about the form of 
the planet spectra, only that it is moving compared to the stationary 
stellar spectrum with specific Doppler shifts. For example, forcing 
the planet spectra, p x (ui) and p y (oj) to be a scaled copy of the stel- 
lar spectra, e(a)s x (u)) and e(a)s y (uj) respectively, would require 
solving three highly non-linear coupled equations, and remove the 
freedom of solving for a wavelength dependent albedo. 



5.1 Fourier Transform 

The Fo urier transform is comp uted using FFTW, a C subroutine 
library l lFrigo & Johnsonl 12005 ffl Before transforming the input 
spectra into Fourier space, the spectra are extended above and be- 
low the wavelength range by adding the value of the continuum, 
and then damped to zero at the ends gradually, using a Hanning 
window which has the form: 



/(A) 



1 / f2n\ 
2{ l - COS { — 



(13) 



for N pixels. FFTW assumes that the spectrum is infinitely periodic, 
and the Hanning window suppresses discontinuities as it loops over 
the wavelength range when transforming into Fourier space. 

The low frequency results of p x (uj) and p y (ui), s x (oj) and 
s y {<jo) need to be dealt with carefully. At these frequencies the 
phase shifts induced by the changes in the planet's velocity be- 
come so small that they are difficult to detect. In essence, the planet 
and star signals are no longer easily separable. Mathematically, this 
problem manifests itself by making the matrix C singular or nearly 
so. As a consequence, the solutions start to be dominated by noise 
or low-frequency systematic errors. By contrast, narrow stellar fea- 
tures that have a presence in both S(A) and P(A) can create measur- 
able amplitudes and phase shifts at high frequencies, and thus they 
contain most of the planet's detectable signal. In order to extract 
this important information, the spectra are effectively put through a 
high-pass filter, and the low frequencies are not solved for. So that 
the result can be inverse Fourier transformed, the low frequency so- 
lutions in Fourier space are replaced with a predicted solution. This 
requires the simple assumption that the planet spectrum is a scaled 
copy of the stellar spectrum. The average of the Fourier amplitudes 
of the input spectra is taken to be an estimate of the Fourier ampli- 
tudes of the stellar spectra. The planet signal will be averaged out, 
due to different Doppler shifts in each exposure. This is multiplied 




. I, . 



Figure 4. Damping spurious signals at low frequencies. Panel A - Exam- 
ple of Fourier amplitudes of input blue mock spectra. Panel B - Example 
of Fourier amplitudes of the solution blue planet (upper subpanel) and star 
(lower subpanel) for mock spectra with e = 0. 1 and low frequencies replaced 
with the artificial signal to damp spurious small frequency signals. For the 
case where the planet-to-star flux is large, the planet and stellar signals ap- 
pear to be the same shape. Panel C - As for Panel B, with e = 0.0001. For the 
smaller e value the moving signal solution is of the order of the statistical 
noise, and hence the planet signal appears to be flat out to large frequencies. 
The central, small amplitude regions in Panels B and C, are due to the need 
for a smooth transition between the solution signal and the artificial signal. 
In this case a Hanning window was used to damp the solution and artificial 
signals to zero where they meet. The frequency range lost in doing this is 
the same for the planet and the star, therefore the resulting planet-to-star 
flux will not be affected. This can be verified with the mock spectra gener- 
ated with particular e values. The Fourier amplitude of the zero frequency 
extends above the scale plotted, but has a finite value. Note that only every 
5th point is plotted here. 



2 See htp://www.fftw.org/ 
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by the predicted amplitude of the planet-to-star flux, and this test 
amplitude is iterated until it equals the solution from the Fourier 
analysis method - the real planet-to-star flux. A null result occurs 
when this iteration does not converge. 

The frequency region of the solution which is artificially re- 
placed is determined for when the matrix C is singular (u> = 0) 
or nearly singular. The gradual transition at the edge of this region 
uses a Hanning window to damp the artificial and solution Fourier 
amplitudes to zero where they meet (see Fig.|4](. This is not the only 
possible method to combine the artificial and real solutions, how- 
ever it must be gradual enough not to cause ringing in the planet 
and star spectra when they are transformed back to real space. Due 
to this transition, sections of frequency space are lost. However, 
as it occurs equally in the stellar and planet spectra (due to the 
matching transition range), the resulting planet-to-star flux value is 
conserved. If the Fourier components are used to compare the stel- 
lar and planet solutions, this artificial replacement is not required, 
however the low frequency result is lost. 

The central frequency cannot just be set to zero to remove the 
low frequency solutions if the inverse Fourier transform is required, 
as it fixes the overall continuum average to also be zero. 

5.2 Measuring the Planet-to-Star Flux 

The planet and the stellar spectra are solved for at each wavelength, 
by taking the inverse Fourier transform of the resulting Fourier 
components p x (oj) and p y {ui), s x (uj) and s y (ui). The spectra 
are retrieved, including the Hanning-windowed continuum-added 
ranges either side of the pixel range corresponding to the solution 
planet and stellar spectra; P(A) and S(A). The planet spectrum so- 
lution is aligned with the stellar spectrum, as it is shifted to the zero 
velocity case. 

The zero frequency value of the spectra is the power in the 
longest wavelength Fourier component, F(0). This corresponds to 
the continuum level, however we cannot explicitly solve for the 
zero shift solution. Therefore, it is artificially constructed before 
the inverse Fourier transform, and is used to iterate to the correct 
solution. The planet spectrum is divided by the stellar spectrum, 
and the mean ratio value taken for the entire range to be the value 
of the planet-to-star flux. This is then tested against the test ampli- 
tude for the artificial low frequency solutions, until the two values 
converge on the solution. 

The planet-to-star flux of the observation is determined with a 
small range of orbital phases, so that the amount of light the planet 
reflects does not vary rapidly. In order to scale this to the maxi- 
mum reflectance for comparison with the Collier Cameron result, 
a planet phase model must be assumed. For the mock spectra, the 
orbital phase is in the range 0.55-0.60. This corresponds to an aver- 
age planet phase of 85 per cent illumination, for a Venus-like phase 
function. The planet phase at each exposure cannot be accounted 
for in solving for the planet spectrum in Fourier space, requiring a 
small region of the orbit to be observed. The actual planet-to-star 
flux e(a) is therefore the model-independent measurement, which 
will be less then the maximum reflectance upper limit. 



5.3 Non-Grey Albedo 

So far in this paper, a grey albedo has been assumed in constructing 
mock spectra. This means that the planet reflects an equal fraction 
of the starlight across the entire wavelength range. In the Collier 
Cameron method, the line weights can be attenuated according to 



a particular atmosphere model to account for a non-grey albedo. 
For exanple, in order to test for different classes of planets, such 
as the Class V roaster or Isolated Class IV models, a wavelength 
dependent planet-to-star flux cold be used for the artifical solution. 
Any increase in the likelihood of the measured planet-to-star flux 
can be used as an indication of the accuracy of the model. 

The Fourier analysis method provides an opportunity to solve 
for the planet signal at each wavelength, and therefore determine 
any slope of the albedo function with wavelength. However, as the 
zero shift solution cannot be solved, the F(0) values are estimated 
based on the assumption that the stellar spectrum can be repre- 
sented by the average input spectrum, and the planet as a scaled 
copy. Therefore, varying the shape of the continuum of the result- 
ing planet spectrum requires modelling possible wavelength depen- 
dent albedos as an iterative parameter to find the best solution. It is 
also possible to divide the input spectrum into short wavelength re- 
gions over which a no n-grey albedo would vary slo wly. This was 
the process adopted by Coll ier Cameron et al.l j 19991) for r Bootis b 
over ranges of 40 nm, and their results suggested a wavelength de- 
pendence of the albedo at around 500 nm. Alternatively, solving for 
the planet's scattered spectrum allows the chance for strong absorp- 
tion lines to be explicitly revealed in high signal-to-noise ratio data. 
Typically these lines are expected to be below the level required 
to detect single absorption lines, but an advantage of the Fourier 
method is that it preserves the details of the planet's spectrum. 



5.4 Monte Carlo Analysis 

Monte Carlo simulations were run for the Fourier analysis method 
with the mock spectra in order to determine the observing strat- 
egy most sensitive to small planet-to-star fluxes, and improvements 
possible with minimal assumptions and typical transit data. 



5.4.1 Signal-to-Noise Ratio 

We constructed spectra with a range of signal-to-noise ratio and 
planet-to-star flux ratios, as per the Monte-Carlo analysis discussed 
in Section 4.1. The results of the Monte Carlo simulations with 
varying signal-to-noise ratios are presented in Fig. [5] correspond- 
ing to the red wavelength region of 560-608 nm and the blue wave- 
length region of 440-515 nm respectively. The planet- to- star flux 
is the oretically predicted to be e < 10 " 4 JCollier Cameron et al.l 
1999); the range of e values is centered on this value. Fig.|2]shows 
the probability distribution and sensitivity level determined from an 
example Monte Carlo simulation, centred on the expected planet- 
to-star flux at maximum reflectance for the system. The resulting 
planet-to-star fluxes are scaled to the maximum reflectance for the 
contour plots presented in Fig. [5] 

The method is more sensitive to smaller planet-to-star fluxes 
in the blue wavelength range than the red wavelength range. This 
is most likely due to the increased number of absorption lines that 
have a larger signal in the Fourier domain, and in the higher fre- 
quency ranges, not as affected by small shifts in the planet spectra. 
For a typical signal-to-noise ratio of around 350, a la detection 
can be made of scattered light with e > 4 x 10 -5 in the blue 
wavelength region, as compared with e > 1 x 10 -4 for the red 
wavelength region of the spectrum. In Fig.[5]the contours of levels 
of detection have a very steep gradient as e decreases to 10~ 5 , sug- 
gesting that very dim planets will require very high signal-to-noise 
ratios to be detectable via this method and ground-based capabil- 
ities. The theoretically predicted planet-to-star flux of e < 10~ 4 
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dCollier Camer on et al .1 19991) is detectable at the 3<j level with rea- 
sonable signal-to-noise values in the blue wavelength region, and 
at the la level in the red wavelength region. Similar to the Col- 
lier Cameron result, matching the la upper lim it of HP 209458b 
measured with the space-based satellite MOST dRowe et"aT]|2008h 
requires the blue wavelength range Subaru HDS spectra with a 
signal-to-noise ratio greater than ~ 900. 



5.4.2 Spectral Resolution 

The spectral resolution of the test input spectra was also varied 
for Monte Carlo simulations with the Fourier analysis method. The 
signal-to-noise ratio was kept constant per pixel, but the oversam- 
pling varied by degrading the spectrum with Gaussian smoothing 
over an increasing radius. The spectral resolution range of R = 
32,000-90,000 corresponds to slit widths of 1.2-0.4" and an over- 
sampling of 4—12 pixels per resolution element. The critical sam- 
pling is 2 pixels per resolution element for Subaru HDS. 

Within this range there is no improvement in the sensitivity of 
the Fourier analysis method. However, increased spectral resolution 
may be beneficial to the data reduction process and correcting for 
instrumental effects and wavelength shifts not incorporated into the 
ideal test spectra. On the other hand, the increase in signal-to-noise 
ratio in echelle data gained by a decrease in spectral resolution may 
yield a measurement of a smaller planet-to-star flux. 



5.4.3 Wavelength Range and Spectrum Length 

The Fourier analysis method was applied to two separate wave- 
length ranges. The two regions of spectra, 560.0-608.0 nm (red 
wavelength range) and 440.0-515.0 nm (blue wavelength range), 
are similar in total wavelength span, however, due to the con- 
stant bin width in logarithmic space, the blue region has many 
more pixels in total. A large number of lines are required to in- 
crease the planet signal abov e the level of the noise dliu et al.l2007l ; 
ICollier Cameron et alj[l999l) . For the Fourier analysis method, the 
blue wavelength range is more sensitive to smaller planet signals 
(see Fig. [5} due to the increased density of sharp absorption lines 
that have a larger signal in the Fourier domain. For the applica- 
tion of the Collier Cameron method to the mock spectra, with ideal 
stellar spectra removal, this is also true. However, for real data, 
where the stellar spectrum is not exactly known, a high density of 
absorption lines in the stellar spectrum is expected to increase the 
difficulty in fitting the continuum and removing the stellar spec- 
trum without affecting the planet signal. For many methods, spec- 
tra with wavelengths towards the red are better suited to detecting 
the planet signal, as absorption lines ar e scattered into regions of 
the spectrum without strong stellar lines Jliu et alj2 007). The con- 
tinuum regions in the red wavelength range also aid the removal of 
the stellar spectrum in the Collier Cameron method, as the template 
spectrum can be more accurately scaled. 

Using a shorter wavelength range is beneficial for decreas- 
ing the effect of assuming a wavelength-independent albedo. The 
Fourier analysis method shows gradual improvement in sensitivity 
to smaller planet-to-star fluxes for spectra covering a larger wave- 
length span, due to the increased number of absorption lines con- 
straining the moving signal. In practice this is limited by the regions 
of bad pixels on a detector, as the Fourier analysis method requires 
continuous spectra. 



5.4.4 Number of Exposures 

For alternative methods of measuring the scattered light contri- 
bution of a non-transiting planetary system such as the Collier 
Cameron method, large numbers of spectra are required to fully 
cover the orbita l phase range of the plan et and determine the 
planet's velocity dCollier Cameron et al]2002T) . This is not required 
for transiting systems, where the velocity of the planet is known for 
the full orbit. As the scattered light signal from the planet is intrin- 
sically dim, taking spectra at maximum planetary phase increases 
the chances of measuring a planet flux. 

Thirty exposures were used to measure the planet-to-star flux 
of the ideal system in testing the Fourier analysis and the Col- 
lier Cameron method. This allows for a range of planet velocities, 
with differing Doppler shifts, but limits the processing time and the 
spread of planet phases. The velocity must be large enough to shift 
the planet signal clear of the corresponding stellar features. Using 
a shorter exposure would yield more images within the possible 
orbital phase range, but reduce the signal in each exposure. 

For up to 50 exposures in a phase range of 0.55-0.60 there is 
an improvement in the sensitivity to smaller planet-to-star fluxes of 
the Fourier analysis method, due to the stronger constraints on the 
moving signal. This results in exposures where the planet velocity 
differs by a couple of m/s in each exposure. 

The number of exposures with sufficiently different veloci- 
ties is limited by the range within which the planet phase does not 
rapidly vary, and the length of the exposure required for an ade- 
quate signal-to-noise ratio. It is important that the separations in 
the planet's velocity from one observation to the next are not pre- 
cisely uniform, as the solutions can blow up when the wavelength 
of the Fourier component is some multiple of the frequency spac- 
ing. The number of spectra used for the Monte Carlo simulations is 
consistent with the number of spectra available for the Subaru HDS 
data of HD 209458, which we discuss in the following Section. 



5.5 Summary 

The Fourier analysis method is generally more sensitive to smaller 
planet-to-star fluxes than the Collier Cameron method for the ideal 
test spectra presented here, due to the ability to extract the planet 
spectrum without having to remove a modelled version of the stel- 
lar spectrum. At the la level, the Fourier analysis and the Col- 
lier Cameron method are similarly sensitive to small planet-to-star 
fluxes in the blue wavelength range. In the red wavelength range, 
the Fourier analysis method is around three times more sensitive 
than the Collier Cameron method at the la level. At the 2 and 3a 
levels, the Fourier analysis method is twice as sensitive as the Col- 
lier Cameron method in the blue wavelength range. The Fourier 
analysis method utilises the known parameters in the transiting 
planet case, and allows a more sensitive measurement from a small 
portion of the orbit that can be observed over a short amount of 
time. This is useful for observations with ground-based telescopes 
for follow up observations to space-based discoveries. 

The Fourier analysis method developed here is best suited to 
studying hot Jupiter planets with orbital periods of a few days, and 
spectra taken of bright host stars in the blue wavelength range, 
just before or after a secondary eclipse. Observations of transit- 
ing hot Jupiter planet's should be planned for when the planet is 
furthest from the observer and therefore appears close to full phase 
and maximum reflectance. The planet should have a velocity large 
enough to shift the scattered spectrum clear of the corresponding 
stellar absorption lines. When the wavelength of the Fourier com- 
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Figure 5. Signal-to-noise ratio for mock spectra with the Fourier analysis method. Contour plot of la, 2a and 3a levels of detection (solid) 
for signal-to-noise ratio values in the range 10 2 -10 3 , and e values in the range 10~ 5 -2 X 10 -4 . Panel A - Results for the red wavelength 
range 560-608 nm. Panel B - Results for the blue wavelength range 440-515 nm. The dots represent the sampling of the parameter space. 
The dot-dashed line shows the signal-to-noise ratio value of 350 for the HD 209458 data. A la level of detection for this data requires a 
planet-to-star flux greater than e = 1 X 10 — 4 in the red wavelength range, and 6 = 4 X 10 — 5 in the blue wavelength range. The dashed line 
shows the upper limit on the planet-to-star flux of HD 209458b set bv lRowe etal] j2008l) of e < 1.6 X 10~ 5 . 



portent is some multiple of a precisely uniform velocity spacing of 
the observations, the matrices become singular and the solutions 
can blow up. Therefore, observations should be planned with some 
intentional irregularity in velocity differences. 

Observations taken over a short velocity range are beneficial to 
reduce the variation due to the unknown planet phase function. The 
relatively small number of exposures required to obtain a measure- 
ment of the planet-to-star flux with the Fourier analysis method is 
therefore beneficial for rapid ground-based follow up to detections 
by space-based surveys. The Fourier analysis method is also suited 
to measuring the planet-to-star flux over short wavelength ranges 
and therefore determining the likelihood of hot Jupiter planets hav- 
ing wavelength dependent albedos in the optical. 



6 SUBARU HDS DATA 

In this Section we use Spectra of HD 209458 and HD 189733 
to test the Fourier analysis method for detecting light scattered 
by a transiting extrasolar planet with ground-based capabilities. 
This data set includes 32 echelle spectra from observations taken 
200 2 October 26 of HP 209458 with Subaru HDS with an 0.8" 
slit dNoguchi et al.|[20o3) and 47 echelle spectra from observations 
taken 2008 August 26 of HD 189733 with Subaru HDS with a 0.4" 
slit. The total wavelength ranges of the CCDs are 554-680 nm for 
the red, and 415-550 nm for the blue when the orders are com- 
bined. 

The HD 209458 data has a spectral resolution of R = 
A/AA = 45,000, and a signal-to-noise ratio per pixel of ~350 
for a 500 secon d exposure. For d etail s of the observatio ns of 
HD 209458b see lWinn et al] d2004» and lNarita etakl d2005l) . The 
HD 189733 data has a spectral resolution of R = A/AA = 90, 000 
and a typical signal-to-noise ratio per pixel of ~300 for a 500 sec- 
ond exposure. During the observations HD 209458b had just left 
a secondary eclipse and had velocities in the range 30-80 km s -1 
and HD 189733b had velocities of 60-130 km s _1 just prior to a 
secondary eclipse. 

The echelle spectra taken were processed using standard IRAF 
procedures, following a correction for the non-linearity of the HDS 
CCDs using the function measured by Tajitsu2010. The one- 



dimensional spectra need to be normalised for variation in the ef- 
ficiency of the CCD between the peak at the centre and the low 
count level close to the edge. An 11-knot cubic-spline was used 
to normalise the large scale variation of the response function of 
the red and blue CCDs, without removing the weak planet signal 
and the absorption lines. The simplest method of combining the or- 
ders of each CCD by averaging the overlapping regions of the nor- 
malised spectra can magnify the noise and was therefore avoided. 
Instead, the object spectra and the response functions were sep- 
arately summed. The combined spectra were then divided by the 
combined response function to remove large scale profiles. This 
preserves the signal in the overlapping regions. Bad pixels in the 
centre of the red CCD and at the edges of the blue CCD limit the 
full wavelength range possible for a continuous spectrum. 



6.1 Correction for Instrumental Variations 

Before the combination of the dispersed orders, the variability 
of the spectra with time due to instrumental effects can be re- 
m oved via a CCD re sponse and wavelength correction as outlined 
in IWinn et al.1 (2004). The mean ratio of two spectra taken at dif- 
ferent times should be equal to unity. This is not the case for the 
Subaru HDS data, as the ratio changes with wavelength bin. The 
spectrum format is affected by instrumental effects such as flexure 
of the spectrograph, changes in the temperature, changes in the in- 
clination of the echelle and cross-dispersion gratings and the flux 
calibration ( Aoki 2002). 

The temperature can change by around 0. 1 degrees per hour, 
causing a 0.14 CCD pixel shift of the spectrum. Changes in the 
collimator mirror or the inclinations of the gratings can cause up 
to 1 pixel difference in the repeatability of the spectrograph set up. 
The largest variation within a night of observations comes from 
the flux calibration, which can vary the profile of the orders of the 
spect ra in different exposures by up to 10 per cent dSuzuki et alj 
120031) . 

The variation between exposures increases with time, such 
that the largest variation is seen b etween the first and last exposure. 
Following the process outlined in Win n et all d2004l) . the pattern of 
the variation of the ratio with wavelength bin can be used to correct 
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adjacent orders to scale the flux calibration to match the first expo- 
sure. The ratio of the n th order of two spectra in the j wavelength 
bin is given by; 

R(Vi) = |4M, (14) 

where R si gnifies the ratio, Si is the first exposure and S2 is a later 
exposure dWinn et alj|2004) . The ratios of the exposures for orders 
adjacent to the one being corrected are boxcar smoothed over j, 
with a width of 100 pixels. The flux corrected n th order of the 
second spectrum S 2 b(A n ) is given by; 

S 2 b(^™j) = [cn+iR'(A„+i,j) + c n _iR'(A„_ij)]S2(A n j), (15) 

where R'(An) signifies the boxc ar smoothed ratio of the n th order 
of the spectra jWinn et al.ll2004h . Linear regression is used to find 
the factors c n +i and c n _i such that the sum of squared residuals is 
a minimum between Si (A„) and the corrected second spectrum. 

A correction can also be made for small variations in wave- 
length. These variations are of the order of a pixel, corresponding 
to a velocity shift of less than 2 kms" 1 , compared to an average 
planet velocity of 70 kms -1 . The n th order of the matched spec- 
trum S2m(A„) is given by; 

dS oh (A) 

S 2 m(A nj ) = c S ab (An,i) + AA ^ ; , (16) 

where the derivative is determined via 3-point Lagrangian interpo- 
lation JWinn et al.ll20oi) . Linear regression is used to find the fac- 
tors Co and AA such that the sum of squared residuals is a minimum 
between Si(A n ) and the matched second spectrum. 

For the Subaru HDS data, orders at wavelengths shorter than 
the location of the bad pixels in the centre of the red CCD were 
corrected, corresponding to a range of 558.83-599.56 nm. The blue 
wavelength range of 445.28-508.11 nm was also corrected before 
normalisation and combination. 



7 RESULTS 

In order to predict the wavelength shift of the scattered spectra with 
respect to the stationary stellar spectra, the velocity of the planet 
was determined from the physical parameters listed in Tableland 
the Julian date at the time of each exposure. During the observa- 
tions HD 209458b had just left a secondary eclipse and had veloci- 
ties toward the observer in the range 30-80 km s" 1 and an average 
planet phase of 0.88. Just prior to a secondary eclipse HD 189733b 
had velocities of 60-130 km s" 1 away from the observer and an 
average planet phase of 0.74. 

The expected level of sensitivity of the Fourier analysis 
method can be determined for the observational parameters and the 
Monte Carlo simulations presented in Section[5] Fig.[5]shows con- 
tours of the level of detection of planet-to-star fluxes in the range 
10~ 5 -2 x 10~ 4 and signal-to-noise ratio values of 10 2 -10 3 for the 
red and the blue wavelength ranges. The dot-dashed line indicates 
the observational signal-to-noise ratio for the Subaru HDS echelle 
sp ectra of HD 20945 8 of 350. The upper limit of e < 1.6x 10" 5 set 
bv lRowe etai] J2008h with the MOST satellite is indicated by the 
dashed line. Upper limits can be derived from non-detections based 
on these Monte Carlo results, however the previous upper limit for 
HD 209458b is below the la level of detection for this data with 
the Fourier analysis method. 

The spectra were analysed in the red (558.83-599.56 nm) 
and blue (445.28-508.11 nm) wavelength ranges separately, as 
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Figure 6. Probability curve for HD 189733 in the red wavelength range. 

P( e truel e spectra)> normalised to an area of unity. The dashed lines show 
the 1(t, 2a and 3cr upper limits on the planet-to-star flux given the results 
of the analysis of the HD 189733 spectra. These correspond to upper limits 
of e < 4.6 X 10 -4 , e < 6.1 X 10 -4 and e < 7.5 X 10~ 4 respectively in 
the wavelength range 558.83-599.56 nm. 

the Fourier analysis method requires continuous spectra. For 
the red and blue wavelength ranges of the HD 209458 data, 
the solutions were e = 7.45 x 10~ 4 and e = -1.50 x 10~ 4 
respectively. A negative planet-to-star flux arising from the noise, 
as seen in the blue wavelength range solution, is not physical 
but may be used to statistically constrain the upper limit on 
the albedo. The red result is larger than expected by theory 
and the upper limit on the p lanet-to-star flux of HD 209458b 
of e < 1.6 X 10~ 5 set by iRowe et al.1 d2008l) . Similarly, for 
the red and blue wavelength ranges of the HD 189733 data, the 
solutions were e = 2.45 x 10 -4 and e = -1.08 x 10~ 4 respectively. 

An upper limit can be set on the planet-to-star flux via proba- 
bility distributions from Monte Carlo simulations and the physical 
constraint that the actual planet-to-star flux must be positive. From 
Bayes' theorem we have: 

P( e truel e spectra) K P( e spectral e true)P( e true)> (17) 

where P(etruei e spectra) ls me probability that the measured 
planet-to-star flux is the physically correct value given the data. 
P( e spectral e true) ls me likelihood that the data will result in the 
measured planet-to-star flux. The prior P(etrue) constrains the 
physical planet-to-star flux to be positive. 

Monte Carlo simulations were run for the simulated data with 
the exact observational parameters of the Subaru HDS data in or- 
der to determine P(espectral e true)> tne likelihood that the data will 
result in the measured planet-to-star flux. For the HD 209458 re- 
sults, 32 spectra were generated with a signal-to-noise ratio of 350, 
a spectral resolution of 45,000 and an orbital phase range of 0.55- 
0.60. For the HD 189733 results, 25 spectra were generated with a 
signal-to-noise ratio of 300, a spectral resolution of 90,000 and an 
orbital phase range of 0.34-0.44. 

The likelihood of the results of the HD 209458 and 
HD 189733 data in the blue wavelength range shows a very small 
probability of having a positive and therefore physical planet-to- 
star flux. This suggests that the data correction did not fully account 
for the time varying signal due to the high density of absorption 
lines in the blue wavelength range, and the results measured with 
the Fourier analysis method are spurious. Similarly, the results for 
the red wavelength range of the HD 209458 spectra are thought 
to be spurious as the likelihood shows a very small probability of 
having a physical planet-to-star flux below the previou sly defined 
upper limit of e < 1.6 x 10~ 5 set bv lRowe et al.1 d2008h . 

Fig. [6] shows the results of the Bayesian analysis for the red 
wavelength range of the HD 189733 data. The measured planet-to- 
star flux was e = 2.45 x 10 -4 . The likelihood from the Monte Carlo 
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HD 209458 system 



HD 189733 system 



Orbital Period 


P orb 


3.52474859 ± 0.00000038 da ys 


2.218573 ± 0.00002 days 


Orbital Inclination 


i 


86.929jl°;° f j , o degrees 


85.79 ± 0.24 degrees 


Mass of Planet b 


M p 


0.64 ± 0.06 Mj 


1.154 ± 0.033 Mj 


Mass of Host Star 


M, 


i 1m +0.066 M 


0.82 ±0.03 M 


Transit Epoch 


To 


2, 452, 826.628521 ± 0.000087 days 


2, 453, 988.80336 ± 0.00024 days 


Semi-Major Axis 


a 


0.04704±0.00048 AU 


0.03099±°;°°°§ 3 au 



Table 2. Physical parameters of the star HD 20 9458 with transiting planet HD 209458 b jKnutson et al. 2007) and of the star 
HD 189733 with transiting planet HD 189733b iBouchv et a l. 2005; Ba kos et alj200d) . 



simulations was fit with a Gaussian distribution and multiplied by 
the prior condition that the planet-to-star flux must be positive. The 
probability that the measured planet-to-star flux is the physically 
correct value given the data, P(etruel e spectra)' i s men normalised 
to have an area of unity. From this curve the la, 2a and 3a upper 
limits on the planet-to-star flux of HD 189733b in the wavelength 
range 558.83-599.56 nm are found to be e < 4.6 x 10~ 4 , e < 
6.1 x 10~ 4 and e < 7.5 x 1CP 4 respectively. 



8 ANALYSIS 

To investigate the level of spurious signals in the data we shuffled 
Doppler shifts with respect to the corresponding exposures. This 
was performed for both stars in both wavelength ranges. A distri- 
bution of planet-to-star fluxes was measured due to spurious time 
variable signals arising out of the instrumental noise. Monte Carlo 
simulations were then used to determine the distribution of false 
planet-to-star fluxes that are not correlated with the Doppler shifts 
of the planet signal. Random velocities were selected in the same 
range as the real data; in the range 30-80 km s _1 for HD 209548 
and 60-130 km s" 1 for HD 189733. This was also done for the 
mock data. 

The distribution of e measure( j values for randomly selected 
Doppler shifts are compared for both wavelength ranges of the 
HD 209458 and HD 189733 Subaru HDS spectra in Fig.0 The dis- 
tributions for the corresponding simulated mock spectra are shown 
in Fig. [8] Although the data corrections for instrumental effects re- 
duce the spread of the spurious signals for the Subaru HDS spectra, 
it is apparent that the time variability of the spectrum has not been 
fully removed. This can be seen by the comparison of the spread 
of the false-planet signals for the mock data, compared to the cor- 
rected Subaru HDS data. 

Similar data corrections based on low order polynomial 
corrections to the CCD response and the cross-correlation of 
spectral lin es have been adop t ed for spectra take n with Sub- 
aru HDS dNarita et alj 120051 ; ISnellen et al] l2008h, the High- 



Resol ution Echelle Spectrograph (HIRES) on Keck l Suzuki et al.l 
l2003h and the High Resol ution Spectrog raph (HRS) on the Hobby- 

Ebe rly Telescope (HET ) dRedfield et al1l2008h . 

lNaritaetal] [2005) found the correction in IWinn et alj l2004h 
to be accurate to approximately the Poisson noise limit for the 
HD 209458 data analysed here, and that the main challenge in 
using ground-based spectra for high precision measurements is 
instrumental variat ion, not the removable telluric contribution. 
ISuzuki et all d2003h calibrate the flux of Keck HIRES spectra to 
within 1 per cent, correcting for vignetting and extraction using ref- 
erence spectra. This suggests that instruments with increased stabil- 
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Figure 7. Probability curves. Histograms from Monte Carlo simulations 
showing the probability of measuring an e measurec j with shuffled planet 
velocities used to obtain the solutions. Panel A - Results for the HD 189733 
Subaru data. Panel B - Results for the HD 209458 Subaru data. Upper sub- 
panels - Results for data in the blue wavelength range. Lower subpanels - 
Results for data in the red wavelength range. The dashed curves show the 
Gaussian distribution with the same mean and standard deviation. 



ity are required for more accurate ground-based follow up charac- 
terisation of extrasolar planets. 

The flux and wavelength calibration correction reduced the 
spread of the spurious signals seen in the distributions in Fig. [7] 
by an order of magnitude. In the mock data the blue wavelength 
ranges have narrower distributions of false planet signals, suggest- 
ing that the increased density of absorption lines in this region is 
more likely to constrain the moving signal to be correlated with the 
planet velocity. This is not seen in the HD 189733 corrected data. 
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Figure 8. Probability curves. Histograms from Monte Carlo simulations 
showing the probability of measuring an « measure d w ' tn shuffled planet 
velocities used to obtain the solutions. Panel A - Results for the HD 189733 
mock data. Panel B - Results for the HD 209458 mock data. Upper sub- 
panels - Results for data in the blue wavelength range. Lower subpanels - 
Results for data in the red wavelength range. The dashed curves show the 
Gaussian distribution with the same mean and standard deviation. 



This may be partly due to the 11-knot cubic-spline chosen to re- 
move the shape of the response of the CCDs being less suitable for 
spectra with a high density of lines. Additionally, in practice, the 
wavelength scale and flux calibration correction may be harder to 
implement for the higher density of spectral lines in the blue wave- 
length range. 

The red wavelength range of the corrected HD 189733 spectra 
is the closest to the mock spectra. It can be seen in Fig.[7]and Fig. [8] 
that the two distributions of false-planets have a similar spread. 
Therefore the time variable signal due to instrumental effects has 
been satisfactorily removed, and the corrected data is comparable 
to the ideal simulated data. From this corrected data the lcr up- 
per limit on the planet-to-star flux in the wavelength range 558.83- 
599.56 nm is e < 4.6 x 1CP 4 . This upper limit is not sensitive 
enough to constrain atmospheric models of hot Jupiter planets or 
rule out a highly reflective upper cloud deck. 

These instrumental variations would also hinder the applica- 
tion of the Collier Cameron method to this Subaru HDS data. It 
is beyond the scope of this paper to apply the Collier Cameron 
method to the non-idealised spectra. Future searches for the light 
scattered by an extrasolar planet will continue to adopt the Collier 
Camer on method, and the method developed bv lCharbonneau et al.l 
(1999). These methods are still required for non-transiting planets. 



For transiting planets, the application of these methods with con- 
strained parameters may yield more stringent upper limits. 

The red wavelength HD 189733 data set showed successful 
implementation of instrumental corrections to Subaru HDS spec- 
tra and the application of the Fourier analysis method. The archive 
data obtained to test the Fourier analysis method did not have suf- 
ficient signal-to-noise ratio to significantly constrain the planet-to- 
star flux of the systems. Improved stability of the instrument while 
observing bright systems with high signal-to-noise ratio will be re- 
quired for follow up characterisation of transiting systems detected 
by space-based missions via the Fourier analysis method. 



9 CONCLUSION 

Measuring planet-to-star flux ratios will continue to guide theories 
on the composition of the upper atmosphere of hot Jupiter planets 
and methods of radiation transfer. This is complementary to the de- 
tections of thermal emission from transiting very hot Jupiter plan- 
ets in the optical and infrared. However, previous upper limits set 
by non-detections of the light scattered by systems such as r Bootis 
( Charb onneau et al.l 199^ ; ICollier Cameron et al.ll999l;lLeigh et all 



l2003ah and HD 75289A JRodler et alj|2008l ; iLeigh et alj |2003bh 
support the need for model independent approaches that exploit 
the benefits of the numerous transiting systems to minimse the un- 
known orbital parameters. It is also advantageous to measure phase 
independent planet-to-star fluxes without assuming a grey albedo. 

In this paper we have introduced a new model-independent 
method for isolating the scattered starlight signal from the host star 
flux in high resolution spectra, that is suited to typical transit data. 
Using mock spectra this Fourier analysis method was shown to 
be more appropriate for typical observations of a well constrained 
transiting system, and therefore more sensitive to smaller planet-to- 
star fluxes, than the currently used Collier Cameron method. Us- 
ing Monte Carlo simulations the sensitivity to dim planets of the 
Fourier analysis method was shown to increase with higher signal- 
to-noise ratio data, and to be better suited to the blue region of the 
visible spectrum due to the increased signal in Fourier space of the 
numerous deep absorption lines. The Fourier analysis method is 
more likely to detect the planet-to-star flux of a dim planet such as 
HD 209458b, with Subaru HDS and ideal high signal-to-noise ratio 
spectra, than the Collier Cameron method. 

The Fourier analysis method for extracting the light scattered 
by transiting hot Jupiter planets from high resolution spectra was 
applied to echelle spectra of HD 209458 and HD 189733. Due to 
instrumental variations that could not be fully corrected for, there 
was no improvement on the measurement of the upper limit of the 
pl anet-to-star flux of HD 209458 compared to the current limit set 
bv lRowe et alj {2008). A la upper limit on the planet-to-star flux 
of HD 189733b was measured in the wavelength range of 558.83- 
599.56 nm of e < 4.5 x 10~ 4 . This result does not constrain atmo- 
spheric models. A deeper measurement of the upper limit of planet- 
to-star flux of this system with ground-based capabilities requires 
data with a higher signal-to-noise ratio, and increased stability of 
the telescope. 

Radial velocity surveys observe stars in a brighter mag- 
nitud e range than transi t ing photometry s urveys teorucki et al.l 
120031: iTinnev et all l200ll ; lUdrv et all [2000). For example Kepler 
is surveying stars of magnitude 14-9, as compared to the Anglo- 
Australian Planet Search (AAPS) on the Anglo-Australian Tele- 
scope (AAT) which targets stars brighter than V< 7.5. Therefore 
measuring the light scattered by the upper atmosphere of plan- 
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ets will most likely require candidates selected via radial velocity 
methods and followed up with transit photometry, to ensure that 
adequate signal-to-noise ratio is obtained. 

The next generation of extremely large telescopes (ELTs) in- 
cluding the European Extremely Large Telescope (E-ELT), the Gi- 
ant Magellen Telescope (GMT) and the Thirty Meter Telescope 
(TMT), will be the largest optical and infrared telescopes ever 
built, dramaticall y increasing th e sensitivity of current ground- 
based telescopes (Carlberg 2005). The reflecting mirror diameters 
range from 25-^2 m. With a diameter of 40 m, the collecting area 
and hence the signal-to-noise will be increased by a factor of ~5 
compared to similar observations with current 8 m ground-based 
telesco pes. This will enable higher sig nal-to-noise ratios to be ob- 
tained dGilmozzi & Spvromiild 12007 ). Therefore, smaller planet- 
to-star fluxes of transiting hot Jupiter planets will be measurable 
with ELTs and the Fourier analysis method. 
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